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Simulating fluid flow in geological formations requires mesh generation, lithology 
mapping to the cells, and computing geometric properties such as normal vec- 
tors and volume of cells. The purpose of this research work is to compute and 
process the geometrical information required for performing numerical simulations 
in geological formations. We present algebraic techniques, named Transfinite In- 
terpolation, for mesh generation. Various transfinite interpolation techniques are 
derived from ID projection op erators. Many geological formations such as the Ut- 
(jTorp arid 



sira formation 



2004: Khattri. He 
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(Maldal and Tappell. 


2004) can be divided into lay- 



ers or blocks based on the geometrical or lithological properties of the layers. We 
present the concept of block structured mesh generation for handling such formations. 
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INTRODUCTION 

Simulation of fluid flow in geological formations, by numerical methods such as Finite 
Elements, Finite Volumes and Finite Differences, requires meshing of the geological for- 
mation into smaller elemen ts called finite volumes or finite elements or cells depending 
on the numerical method (IKhattri et al" . 20061 : Ewing and Heinemann . 1984 : Khattri . 



20061 200,4 Khattri and Aavatsmarkl . l2006| ). These elements in three dimensions can be 



hexahedra, tetrahedra, prism and pyramid. In this paper, we focus only on hexahedral 
mesh generation. It is desirable that the part of the geological formation where solution 
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shows nonlinear changes shou ld be refined ( Khattri and Aavatsmark . 20061 : Khattri . 2005; 



Khattri and Fladmark , 20061 ). Such a sol ution behaviour can occur due to litholog i cal or 



geometrical properties of the formations ( Khattri and Aavatsmark , 20061 : Khattri , 20051: 
Khattri and FladmarkL l2006h . 



Many geological formations and reservoirs of interest can be divided into layers based 
on the geological characteristics such as faults and pinchouts or the lithological properties 



on tne geological cnaractenstics sucn as iauits ana pmcnouts or tne ntnoiogicai properties 
such as shale and s andstone. For example, th e Utsira formation dTorp and Galel . 2004; 



Khattri et al. . 20061 ) and the Sn0hvit gas field ( Maldal and Tappel . 12004 ) . Each of these 



layers can be meshed into hexahedrals by the algebraic techniques independent of the 
other layers. In this way grid distribution and quality of mesh can be improved and 
controlled in each of the layers separately. This technique is called the multilayer or 
the multiblock approach. The concept of multiblock mesh generation is very useful for 
handling layered formations. Some of the advantages of this approach are 

1. Many geological formations can be realized by this concept. 

2. It makes parallelization of a single phase problem straight forward. The multi- 
block/multilayer approach used as a domain decomposition concept allows the direct 
parallelization of both grid generation and flow codes on massively parallel systems. 

3. Grid density, distribution and quality can be controlled easily. It is desirable that in 
the areas of expected great nonlinear changes of solutions (around wells and material 
discontinuity) mesh should be refined. 

4. Controllability over the simulation. For example, the implementation of lithology 
and local optimization of mesh quality. 

5. Though at the global level multilayer grids are unstructured in nature. Still at local 
level mesh can be expressed by logical numbering. Optimization of the quality of 
structured grids is easier. Instead of performing global mesh optimization, mesh can 
be optimized around critical locations such as wells. Structured grids can easily be 
made orthogonal at the boundaries and also almost orthogonal within the solution 
domain thus facilitating implementation of boundary conditions and also increase 
numerical accuracy. Discretization of partial differential equations on structured 
meshes is easier than on unstructured meshes. 

6. A structured grid produces a structured matrix and thus makes it easier to use 
sophisticated linear solvers. 

Now let us discuss about algebraic method of grid generation. 

Algebraic Method of Mesh Generation 

In the algebraic method of grid generation, we seek an algebraic mapping from a cube 
in computational or reference space to a physical space with the corresponding boundary 
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Figure 1: Mapping a unit cube onto a physical domain. 



surfaces ( Knupp and Steinberg . 1994 : Knuppl . 1991 ). Transfmite interpolation (TFI) is 
such an algebraic mapping. TFI is also referred to as multivariate interpolation or Coons 
Patch. Figure ^ shows a mapping from a unit cube in the reference space onto a physical 
domain. Let the reference or computational space be defined by £, r] and k coordinates, 
and the physical space be defined by x, y and z coordinates. Suppose there exists a 
transformation or mapping, r = r(£, t],k), which maps the unit cube onto the interior 
of the physical domain, and this mapping maps the boundary surfaces of the cube to 
the corresponding boundary surfaces of the physical domain. Thus, rj — 1 surface of the 
cube is mapped to the r(l, tj,k) boundary surface of the physical domain. Transfmite 
interpolation is the boolean sum of univariate interpolations in each of the computational 
coordinates. Univariate interpolations are also referred to as one dimensional projection 
operators or projectors. Boolean sum of the projection operators are defined below. 
A univariate interpolation is an operator that vary only in one dimension or roughly 
speaking it is a function of only one reference coordinate. A univariate interpolation can 
be linear, quadratic and cubic. Any univariate interpolation can be applied in a coordinate 
direction. Generally a higher order interpolation operator is desired in flow direction. TFI 
is composed of ID projection operators, let us first define some one dimensional projection 
operators. 



One Dimensional Projection Operators 

A ID projection operator or projector can be defined in many ways depending upon the 
available information. For example, a linear projector can be formed from two surfaces; 
a Hermite projector can be formed from two surfaces and directional derivatives at these 
surfaces; a Lagrangian projector can be defined from two boundary surfaces and internal 
surfaces. 

Let the reference space be defined by £, r] and k coordinates (£ G [0,1], rj G [0,1] 
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and k G [0, 1]). Suppose there exists a transformation r(£,rj, k) from a unit cube in the 
reference space onto a physical domain. That is r: k i — ► k. Let the physical space 
be defined by six boundary surfaces. A £ surface in the physical space is a surface on 
which value of £ is constant. Thus, two £ boundary surfaces are r(0, 77, k) and r(l,7/, ft). 
Similarly, two 77 and k boundary surfaces are given as r(£, 0, k), r(£,l,/c) and r(£,77, 0), 
r(£, 77, 1), respectively. From these six boundary surfaces, the following ID projection 
operators are defined 

P e = f (l-Or(0,77,K) + £r(l,77, K ) , (1) 
P r? H f (l-r7)r(e,0,K)+r7r(e,l,K) , (2) 
P K d ^ f (l-/ t )r(£,r/,0) + K r(e,r / ,l) . (3) 

The projectors Pg, P^ and P K are ID projection operators and they are functions of the 
coordinates (£,77, ft). The projection operators defined by equations (0), (0) and Q are 
linear in £, 17 and ft coordinates. It can be notice that the operators are defined from two 
surfaces of a particular kind. For example, P^ is defined from two £ boundary surfaces in 
the physical space r(0, i], ft) and r(l, i], ft). 



If in addition to the boundary surfaces we also know the internal surfaces of a domain 
then a projection operator can also be defined from more than two surfaces of a kind. 
For example, if there are n + 1 surfaces of £ type (n — 1 internal curves and 2 boundary 
surfaces) then Pg projection operator can be defined as 



def 



(4) 



j=0 



fjBerrut and Trefethenl . l2004t iHighaml . I2004T ) . Here, j — and j = n are the boundary 
surfaces while j = 1, . . . , n— 1 are the internal surfaces, and (3j is the Lagrangian weighting 
factor. The Lagrangian weighting factor is given as follows 



ft<0= n (M) 



(5) 



It can be notice that the weighting factor /%(£) is an order n polynomial having zeros at 
all of the surfaces except the jth surface. The Lagrangian weighting factor satisfies the 
following 

Now let us express the Lagrangian projection operator in another form. The numerator 
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in the Lagrange weighting factor (jHJ) can be written as 

[(£ - GO K - 6) • • • (£ - 60] 



(7) 



(see 



2mM 



3errnt and Trefethen , 20041 ) . Let us define the barycentric weights as ( Berrut and Trefethen , 

1 



Usi ng equations (17|) and (151). th e Lagrangian projection operator 
as ( Berrut and Trefethen . 200j 



def 



j=0 



(8) 

can also be written 
(9) 



In the grid generation litera ture, the equation (El) is used but the new form (jHJ) is compu- 
tationally more efficient (cf. iBerrut and Trefethenl 12004^ . Similarly, if in addition to the 
boundary surfaces we are also given the derivatives (direction vectors) on these boundary 
surfaces then we can define the Hermite interpolation operators. For example, if we are 
given two £ surfaces: r(0, rj, k) and r(l, r], k), and let the direction vectors on these sur- 
faces be r'(0, rj, k) and r'(0, rj, k), respectively. Then, the ID Hermite projection operator 
can be defined as 



P 5 = f (2£ 3 



3 e + 1) r(0, r,, k) + (-2 £ 3 + 3 £ 2 ) r(l, r], «) 



Hermite projectors are easy to implement and are powerful tools for grid generation. Grid 
lines can be made orthogonal by the proper choice of direction vectors. This may help in 
accurate modelling of boundary conditions. 



Figure 12 shows a physical domain containing three £, three rj and two k surfaces. Since 
the domain contains three £ surfaces, three rj surfaces and two k surfaces thus we can 
define a Lagrangian operator, a Lagrangian operator and a linear P K operator. 
Figure El shows another physical domain with two £, three rj (r(£, 0, k), r(£, r/i,K) and 
r(£, 1, k)) and two k surfaces. For this domain, a linear P^, a Lagrangian P^ and a linear 
P K operators can be defined. For this domain, the Lagrangian P^ operator is given as 



p v = n 



T]-0 
where Q is given as, 



^ r(e, 0, «) + ( J*- ) r(£, m , K )+(J^-) r( e, 1, «) 



J} - Vi. 



77 — 1 



'IF 
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Figure 2: A 3D physical domain con- 
taining 3 £, 3 i] and 2 k surfaces. 




77 — 1 



77 = 771 



Figure 3: A 3D physical domain con- 
taining 2 £, 3 17 and 2 k surfaces. 



and c^o, u>i, u 2 and co> 3 are given as, 
1 



u 



1 



-771) (1 — 771, 



UJ 2 ~- 



(-1) (m - 1) 



(12) 



Now let us study two important and useful properties of projection operators. These 
properties are called tensor product and boolean sum of projection operators. 

Properties of Projection Operators 

This section presents two important properties of projection operators. 

0.1 Tensor Product 

Tensor product of the projection operators and P^ is defined as follows 



p fo , d =^ f p f o p„ = (1 - [p„k =0 + i [p„k =1 



(13) 



Here, P^ is assumed to be linear projection operator as defined by the equation (JTJ). It 
is clear from equation (JTHJ) that P^ is a projection operator. That is P^og = P$- If P§ is 
Lagrangian projection operator then the tensor product is defined as 



(14) 



3=0 



Tensor product of two projection operators is also a projection operator (Pg 0f? is a pro- 
jection operator). Since tensor product is also a projection operator, it is commutative 
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in nature. That is = P^- Similarly tensor products can be defined for an arbi- 
trary number of projection operators. For example, the tensor product of three projection 
operators is defined as follows 

P^o K = f P C O (P„ O P K ) = (1 - [Pr pK h=0 + £ [Pr P nh=l ■ (15) 



In the above equation, the projection operator P^ is linear. 



0.2 Boolean Sum 

Boolean sum of two projection operators is a also a projection operator and it is defined 
as follows 

P?0r? = P$ ffi P»7 = + P v — Pgor? • (16) 

Here, P^ is the tensor product of the P^ and P v projection operators. Boolean sum 
is commutative in nature. That is P^ ffi P,, = P,, ffi P ? . Since boolean sum is also a 
projection operator thus it follows the projection property. That is P^©£ = P$. Similarly, 
the boolean sum can also be defined for an arbitrary number of projection operators. The 
boolean sum of three projectors is defined by using the fact that boolean sum and tensor 
product of two projection operators are also projection operators. Thus, the boolean sum 
of P^, P^ and P K operators is given as 

P^©»?©k = P? © Pr? © Pk , 

= P € ©(P,©P K ) , 

= P^ ffi (P-q ffi P K — P^ok) j 

= P^ © Pr/ + P( © Pk P^ ffi PryoK , 

P^ "I - P-r) P^or; ffi P^ ffi Pk P^ok P^ P^ok ffi P^oryoK • 



Thus, 



P^ ffi Pr; ffi Pr P^o»? P^ok Pjyo/i ffi P^ojyoK • (1^) 



Here, P^ denotes the tensor product of Pg and P^ projection operators and P^ok 
denotes the tensor product of P^, P v and P K projection operators. 



Transfinite Interpolation 

Boolean sum of projection operators is the basis for Transfinite Interpolation. TFI are 
extensible used for algebraic grid generation. Since, ID projection operators comes in 
many flavours such as the Lagrangian and the Hermite thus TFI can be defined by many 
different expressions depending upon which ID projection operators are used. Linear 
Transfinite Interpolation creates a grid in 3D using surfaces that define the boundaries. 
Quality of the generated grid strongly depends on the parametrizations of the boundary 
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curves. In its simplest form this mapping blends two given surfaces to create a grid in 
the region bounded by the surfaces or curves. Linear Transfmite Interpolation mapping 
is defined from six surfaces. Transfmite Interpolation mapping will only give a reasonable 
grid if the surfaces that define the boundary match at the edges, and the surfaces are 
parametrized in the same direction otherwise grid lines could cross each other. We are 
using the equation (|17|) for mesh generation. Thus, the position vector in the physical 
space is given as 

r(£, n, k) = P ee „© K = P 5 © P„ © P K . (18) 

Let the geological formation be defined by the six boundary surfaces r(0, rj, k), r(l,n, k), 
r(£, 0, k), r(£,l,/c), r(£, ?y, 0) and r(£, rj, 1). Thus, from these six boundary surfaces the 
linear projection operators can be defined. Let us divide the reference unit cube into nx 
subdivisions in the £ coordinate direction, ny subdivisions in the rj coordinate directions, 
and nz subdivisions in the k coordinate direction. Thus for this mesh 

Number of nodes =nx x ny x nz, 
Number of cells =(nx + 1) x (ny + 1) x (nz + 1), 
Number of surfaces =nx x ny x (nz + 1) + nx x nz x (ny + 1) + ny x nz x (nx + 1). 

A simple routine for generating mesh in the geological formation is given as 



Algorithm 1: Grid generation in a block or layer. 

1: for (ix = 0; ix < nx + 1; ix ++ ){ / / Moving in the £ direction 

2: for (iy = 0; iy < ny + 1; iy ++ ){ / / Moving in the n direction 

3: for (iz = 0; iz < nz + 1; iz ++ ){ / / Moving in the k direction 

4: i := ix + (nx + 1) x iy + (ny + 1) x iz; // Node number 

5: £i := ix/nx; n\ := iy/nx; K\ := iz/nx; / / Gridding of Unit Cube 

6: r(ix, iy, iz) := [P^® K ] ?=Cl v=r]1 K=Kl // Position in the Physical Space 

7: } 

8: } 

9: } 



Computing Geometric Properties 

Let us consider the steady state pressure equation of a single phase flowing in a porous 
medium 

-div(Kgradp) = / . (19) 

In porous media flow, the unknown function p = p(x, y) represents the pressure of a 
single phase, K is the permeability or hydraulic conductivity of the porous medium, and 
the velocity u of the phase is given by the Darcy law as: u = — l^gradp. For solving 
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Figure 4: Division of a Figure 5: Volume of the Figure 6: Normal vector 
hexahedra. tetrahedra 1245. on the surface 1234. 

partial differential equations (PDEs) in geological formations by numerical methods such 
as Finite Volumes, the domain is divided into smaller elements. The process of dividing 
geological formations into smaller elements is referred to as meshing of the domains or 
geological formations, and the elements are called finite volumes or cells. Integrating 
equation ()19j) over one of the finite volumes with volume Vol and boundary <9Vol, and 
using the Gauss divergence theorem leads to 

- / KVp-h = [ f , (20) 
J a vol J vol 

where n is the outward unit normal on the boundary <9Vol of the finite volume Vol. Let 
us assume that finite volumes are hexahedras. Boundary of these finite volumes consists 
of six surfaces <9Volj. The above equation can be written as 

-J2 [ KWp-h= [ f , (21) 

i=l JdVoli -'Vol 

the term — J 9Vol . K Vp ■ n is referred to as the flux or the Darcy flux through the surface 
<9Volj. The term J Vol / can be approximated as value of the function / at the center of 
the hexahedra times the volume of the hexahedra. Thus, converting a partial differential 
equation into an algebraic equation requires volume of hexahedra and normal vectors on 
the surfaces of the hexahedra. Now, let us present a method for computing volume of the 
hexahedra. 

Figure E] shows a hexahedra 12345678. Let the position vector of the vertix i be r, 
with i — 1, ... ,8. This hexahedra can be divided into two prisms 124568 and 134578. 
Each of these prisms can divided into three tetrahedras. The Figure 0] shows the division 
of the prisms 124568 into three tetrahedras 1245, 2456 and 4568. Thus, a hexahedra 
can divided into six tetrahedras, and the volume of the hexahedra can be computed by 
summing the volume of the six tetrahedras. Figure El presents the tetrahedra 1245. The 
vectors V 1; V2 and V 3 are meeting at the vertix 1 of the tetrahedra. The vectors V 1; 
V 2 and V 3 are given as V[ = r 2 — r 1; V 2 = r 4 — r! and V 3 = r 5 — r 1; respectively. The 
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Figure 7: A multiblock grid in a geological formation. 



volume of the tetrahedra 1245 is given as 

Vol 1245 = ^ | Vi • (V 2 x V 3 )| . (22) 
o 

Now, we are going to see two techniques for computing normal vectors on the surface of 
hexahedra. 

For the surface 2468, see Figure 13 The diagonal vectors V 28 and V 46 of the quadri- 
lateral surface 2486 of the hexahedra are given as V 2 § = r§ — r 2 and V46 = r$ — r^. 
The normal vector on the quadrilateral surface is given as the cross product of these two 
diagonal vectors. That is n = V 28 x V 46 . 

The position vector of a point in the physical space (geological formation) is given 
by the expression (|18|) . and this expression is a function of the coordinates £, 77 and k. 
Differentiating this expression with respect to a particular coordinate will give us a vector 
pointing in that coordinate direction. This vector is called the covariant vector. Figure |H1 
presents two covariant vectors and r K . Differentiating the expression (fTHj) with respect 
to 77 results 

^ 9Prj dP £or] dP-qoK ^ dP £ r)OK (23) 

v drj drj drj drj 

Since, and P K are not functions of 77 so their differentiation with respect to 77 will 
vanish. Similarly, the covariant vector r K can be determined. Cross product of these two 
covariant vectors will provide the normal vector on the surface. 



Example 

The geological formation is shown in figure (J7|) is divided into nine layers based on the 
medium property. Four of these nine layers are highly permeable thus these layers are 
densly meshed, as shown in the figure [7| 



10 



ACKNOWLEDGEMENTS 



We thank Ivar Aavatsmark for providing useful comments, and Many L. Buddie and 
David R. Wood for correcting the manuscript. 

References 

Berrut, J.-P. and Trefethen, L. N., 2004, Barycentric Lagrange interpolation: SIAM Rev., 
v. 46, no. 3, p. 501-517 (electronic). 

Ewing, R. E. and Heinemann, R. F., 1984, Mixed finite element approximation of phase 
velocities in compositional reservoir simulation: Computer Methods in Applied Me- 
chanics and Engineering, v. 47, no. 1-2, p. 161-175. 

Higham, N. J., 2004, The numerical stability of barycentric Lagrange interpolation: IMA 
J. Numer. Anal., v. 24, no. 4, p. 547-556. 

Khattri, S. and Aavatsmark, I., 2006, Numerical convergence on adaptive grids for control 
volume methods, Submitted in a Journal. 

Khattri, S. K., 2005, Numerical Analysis of an Adaptive Finite Volume Method for Single 
Phase Flow in Highly Heterogenous Medium, Submitted in a Journal. 

Khattri, S. K., 2006, Analyzing Finite Volume for Single Phase Flow in Porous Media, 
Accepted in the Journal of Porous Media. 

Khattri, S. K. and Fladmark, G. E., 2006, Which are Better Conditioned Meshes Adap- 
tive, Uniform, Locally Refined or Localised: In the 6th International Conference on 
Computational Science, The University of Reading, UK. 

Khattri, S. K., Hellevang, H., Fladmark, G. E. and Kvamme, B., 2006, Simulation of 
long-term fate of C0 2 in the sand of Utsira, Submitted in the Journal. 

Knupp, P. and Steinberg, S., 1994, Fundamentals of grid generation: CRC Press, Boca 
Raton, FL, with 1 IBM-PC floppy disk (3.5 inch; HD). 

Knupp, P. M., 1991, Intrinsic algebraic grid generation, in Mathematical aspects of nu- 
merical grid generation, SIAM, Philadelphia, PA, v. 8 of Frontiers Appl. Math., p. 
75-97. 

Maldal, T. and Tappel, I. M., 2004, C0 2 underground storage for Sn0hvit gas field devel- 
opment: Energy, v. 29, no. 9-10, p. 1403-1411. 

Torp, T. A. and Gale, J., 2004, Demonstrating storage of C0 2 in geological reservoirs: 
The Sleipner and SACS projects: Energy, v. 29, no. 9-10, p. 1361-1369. 



11 



